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Abstract 

We show that deterministic systems with strong nonhnearities seem to be more 
appropriate to model sandpiles than stochastic systems or deterministic systems 
in which discontinuities are the only nonlinearity. In particular, we are able to 
reproduce the breakdown of Self-Organized Criticality found in two well known 
experiments, that is, a centrally fueled pile [Held et al. Phys. Rev. Lett. 65 (1990) 
1120] and sand in a rotating tray [Bretz et al Phys. Rev. Lett. 69 (1992) 2431]. 
By varying the parameters of the model we recover Self-Organized Criticality, in 
agreement with other experiments and other models. We show that chaos plays a 
fundamental role in the dynamics of the system. 
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In 1987 Bak, Tang and Wiesenfeld showed that certain open, dissipative, spa- 
tially extended systems spontaneously achieve a critical state characterized by 
power-law distribution of event sizes(|l|). They denoted this phenomenon Self- 
Organized Criticality (SOC) and illustrated it using a simple two-dimensional 
cellular automaton model for sandpiles. An enormous amount of research fol- 
lowed that work. Some natural phenomena that have been connect to SOC 
are earthquakes, evolution, interface dynamics, solar flares, vortices in super- 
conductors, charge density waves among others. For a review see(j3;0). 

The first observation of SOC in granular material was made by Frette et. 
ai(|3). They observed that elongated grains of rice dropped between two par- 
allel plates with a narrow separation between them (that is, in a quasi-one- 
dimensional geometry) present avalanche distribution consistent with SOC. 
Using more symmetric rice a stretched exponential was seen. 
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For systems in which the flow of grains is on a two dimensional surface different 
behaviors have been reported. To our knowledge, the flrst of such experiments 
was performed by Jaeger et al.Q). They showed that avalanches of sand in 
a rotating drum do not present power-law distribution of event sizes, instead 
they found system-spanning avalanches with a narrowly peaked distribution. 
Similar results were found by them for a avalanches of sand in a box with 
an open side(jl). Using a different set up. Held et ai.(0) studied avalanche 
distribution by fueling a canonical pile on its apex. They observed that if 
the base of the pile is smaller than a given size, then the avalanche distribu- 
tion is consistent with SOC. If the base of the pile is increased beyond that 
critical value, another regime dominated by big, nearly periodic avalanches 
is observed. Bretz et al. performed experiments with granular material in a 
rotating tray(0). They found big spanning avalanches that reset the system, 
but between those large avalanches a sequence of small avalanches was found 
which presents power-law distribution. 

New experiments on granular materials with flow in two dimensions have 
shown avalanche distribution fully consistent with SOC. The main difference 
of these experiments with respect to the previous ones was a change in the 
grain type. Costello et al. studied the avalanche distribution in a canonical 
bead pile(|^ whereas Aegerter et al. studied a rice pile in a box where grains 
are dropped at one of its edge(j^. In both experiments power-law distribution 
of event size and flnite size scaling data collapse was seen. 

Recently, two simple nonlinear two-dimensional deterministic models governed 
by coupled maps lattice (CML) have been introduced by us which repro- 
duce the avalanche distribution observed in the experiments with granular 
materials (fioh . We modeled the experiments in which a canonical pile is fueled 
at its apex and the experiment in which the granular material is in a rotat- 
ing tray. We observed that in both models there is a region of the parameter 
space where SOC occurs, and another region where SOC breaks down. The 
way SOC breaks is different from one model to the other and it is as in the 
experiments. 

Here we expand that work to show that the universality class of those systems 
is not unique. That is, in each model it changes with the parameter values. We 
also study the Liapunov exponent of the models and show that at the onset of 
exponential instability the distribution of event sizes crosses from a power-law 
one to one dominated by big avalanches. Strong chaos is therefore a deciding 
factor for SOC to be present or not, and is closely related to the amount of 
friction between the grains. 

For the sake of completeness we review the algorithms for evolving our systems(0). 
We flrst describe the one that governs the centrally fueled pile, called in (fl^ 
"local dropping" (LD): (1) Start the system by assigning random initial val- 
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ues for the variables Xij, where Xij is the local slope of the pile, so the they 
are below a chosen, fixed, threshold Xth- (2) Choose a nearly central site of 
the lattice and update it slope according to Xij = Xth- (3) Check the slope 
in each element. If an element i,j has Xij > Xth, update Xij according to 
x'i j — (t>{xi^j — Xth), where is a given nonlinear function that has two pa- 
rameters a and d. Increase the slope in all its nearest neighboring element 
according to x'^^ = Xnn + where Ax = Xij — x'^j and nn denotes nearest 

neighbors. (4) If x'^j < Xth for all the elements, go to step (2) (the event, or 
avalanche, has finished). Otherwise, go to step (3) (the event is still evolving). 
Without losing generality, we can take Xth = 1. In our simulation in step (2) 
we have chosen the site with i — j — L/2, where L is the lattice size. The 
nonlinear function we use is 

'^(^^ = |o, otherwise. 



The parameter d is in the interval (0; 1] and is associated with the minimum 
drop in energy after an event involving one single element. The parameter a is 
greater or equal to zero and is associated with the amount of dynamic friction 
between the grains. That is, the smaller the a, the larger the friction and the 
smaller the change in the slope of the pile. 

The algorithm that governs the model used for rotating tray experiment, that 
is, the "uniform driving case" (UD), is similar to the one given above, with 
the exception of step (2). That step is replaced by: (2) Find the element in 
the lattice that has the largest x denoted here by x^ax- Then update all the 
lattice elements according to Xij — > Xi^j + Xth — x^ax- 

As one evolves the system one can study several quantities, such as avalanche 
size (which is the number of site updates in a avalanche), time duration, radius 
of gyration, etc. Here we limit our study to the avalanche size distribution. We 
find that in both models, when d = 1 the only nonlinearity in the system is a 
discontinuity between Xth and 0(0) and we observed that in this case SOC is 
always present. Note that when d = 1 the UD case reduces to the conservative 
OFC model, which is known to present SOC. When a < 1 SOC is found 
for any value of d, whereas when a > 1 and d < 1 scaling consistent with 
SOC happens only for L below a given value. For L larger than that critical 
size, big avalanches that belong to a different distribution arc observed. This is 
consistent with what is found in experiments with granular materials. In Fig.l 
we show the probability distribution of avalanche sizes P{s) as a function of it 
size s for a regime with small d and a > 1 (d = 0.1 and a — 1.2). In (a) we show 
the case of local dropping for L ranging from 8 to 256. In (b) we show similar 
distribution for the uniform driving case for the same system sizes. One sees 
that for small systems sizes (L < 32 in the LD case and L < 16 in the UD case) 
a scaling regime consistent with SOC is seen. We have found finite size scaling 
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Fig. 1. Distribution of event sizes for a = 1.2, d = 0.1, L = 8, 16, 32, 64, 128, 256 for 
the (a) LD and (b) UD cases. 



data collapse for both cases and for the LD case it is reported in ()10l). This 
scaling regime reproduces what is observed experimentally in the pile fueled 
at its apex(0). As the system size is increased a new regime is seen where the 
SOC scaling is lost by the appearance of frequent large avalanches, consistent 
with experiments. Note that for the regime of broken SOC the distribution of 
the small events do not change as the system size increases in the LD case, but 
the the UD case the small avalanches become less frequent as the system size 
increases. This could explain why in some experiments, such as ((5|) modeled 
by the UD case only large avalanches was seen. We have also found that the 
big avalanches in the UD case involves all the elements of the system, and this 
does not necessarily occurs in the LD case. 

We now turn our attention to the regime when a < 1, where the system 
presents SOC for any value of d and L. We plot in Fig.2(a) the avalanche 
distribution for the LD case for two different values of d, that is, d = 0.1 
(solid line) and d = 1 (dashed line) with a = 0.9. One sees that -P(s) ~ s~'^, 
characteristic of SOC systems but the value of r is different from one case to 
the other. We have varied the system size and used finite size scaling ansatz 
by plotting P{s)s'^ = s/L^ for L = 64, 128,256. The results are displayed in 
Fig. 2(b). The lower curves are for d = 1 and the upper ones for d = 0.1. We 
observe that finite size scaling is reasonably well obeyed except at the region 
where the curves bend for d = 0.1. We have used r = 1.26 and D = 2.6 for 
d = 0.1 and r = 1.33 and D = 2.8 for d = 1. We see that the critical exponents 
T and D vary as we vary the parameters. Therefore, the universality class of 
this system is not unique. In Fig. 2(c) we show similar results for the UD case. 
The solid line is for d = 0.1 and the dashed one for d = 1. In Fig. 2(d) we show 
the data collapse using the finite size scaling ansatz for L = 64, 128, 256. The 
curves with a smoother bending are for d = 1. We have used r = 1.33 and 
D = 2.7 for d = 0.1 and r = 1.27 and D = 3.0 for d = 1. Again, we see a 
change of universality class when the parameters of the system are varied. 

We have studied the largest Liapunov exponent (LLE) of the system to verify 
the presence or absence of chaos. If the LLE is positive, the system is said to 
be chaotic. To study the LLE we use the algorithm developed by Benettin et 



ai. ljlll ) and the results are displayed in Fig.3. In (a) and (b) we show the LLE as 



a function of a for the LD and UD cases, respectively. There we have d = 0.1, 
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Fig. 2. Distribution of event sizes for (a) a = 0.9, L = 256 with d = 0.1 (solid) and 
d = 1 (dashed) for the LD case, (b) Finite-size scahng data collapse for a = 0.9, 
L = 64, 128, 256 with d = 0.1 (upper set of curves) and d = 1 (lower set of curves) 
for the LD case, (c) and (d) Same as (a) and (b), respectively, except that now it 
applies for the UD case. 

L = 32 (dashed) and L = 64 (solid). The figures show that the LLE is positive 
for any value of a and that at a = 1 a transition happens. If no discontinuities 
existed between the Xth and 0(0) two trajectories with slightly different initial 
conditions would tend to converge when a < 1 (see Ref.((lO; E3)) resulting 
in a negative LLE. However, because of that discontinuity two competing 
tendencies exist: one for convergence and another for separation. However, 
the one for separation is stronger and this results in a positive LLE, as Fig. 3 
shows. This is a weak form of chaos. When a > 1 there are no competing 
tendencies, and the trajectories always separate, and this will obviously result 
in a positive LLE. This is a strong form of chaos. The transition a < 1 (where 
SOC is found) to the one in which a > 1 (where SOC is broken for large L) is 
markedly different in the two models. 

How the LLE vary with the system size is shown in Fig. 3(c) for the LD case 
and Fig. 3(d) for the UD case. The dashed line refers to a = 1.5 and the solid 
one to a = 0.5. There, d = 0.1. One sees that for small L in both curve the 
LLE tends to decreases. However, as L grows the case in which a < 1 (where 
SOC is present for any L) the LLE continues to get smaller and goes to zero 
(as a power-law function for the LD case), whereas when a > 1 the curves 
change their behavior as L gets bigger (which is consistent with the fact that 
SOC breaks for a given value of L when a > 1). That is, as the system size 
increases the chaotic behavior of the system decreases rapidly in the SOC 
regime, whereas this does not occur in the non SOC regime (in this regime 
the LLE tends to increase or stays approximately constant). 

We summarize our results as follows: the introduction of strong nonlinearities 
in SOC systems causes the appearance of a new behavior as the parameters 
are changed. This is not seen in SOC systems in which discontinuities are the 
only nonlinearity. That is, one is able to see the breakdown of SOC as the 
system size increases, which is seen in experiments. Chaos play an important 
role, since it is only in the regime of strong chaos that the breakdown of SOC 
occurs. 



5 



Fig. 3. The largest Liapunov exponent as a function of a for d = 0.1, L = 32 
(dashed) and L = 64 (sohd) for (a) LD and (b) UD cases. The largest Liapunov 
exponent as a function of L for a = 0.5 (solid) and a = 1.5 (dashed) for (c) UD and 
(d) UD cases. 
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